Scaling and universality of multipartite entanglement at criticality 
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Using the geometric entanglement measure, we study the scaling of multipartite entanglement 
in several ID models at criticality, specifically the linear harmonic chain and the XY spin chain 
encompassing both the Ising and XX critical models. Our results provide convincing evidence that 
ID models at criticality exhibit a universal logarithmic scaling behavior ~ log 2 ^ in the multipartite 
entanglement per region for a partition of the system into regions of size £, where c is the central 
charge of the corresponding universality class in conformal field theory. 



The study of the connection between quantum entan- 
glement and the properties of spatially extended many- 
body systems such as spin chains [3, 0, Q an( i har- 
monic chains 0, [H , has recently attracted considerable 
attention. This connection is especially relevant for sys- 
tems near, or at, certain quantum phase transitions, 
where well-studied features of criticality such as scale- 
free behavior and universality are in fact manifestations 
of the entanglement properties of the underlying quan- 
tum state. 

Universal entanglement signatures of criticality have 
by now been well established in the case of pure, bipar- 
tite entanglement. Specifically, for ID critical systems, 
the entanglement entropy of a region of size £ and its com- 
plement is seen to follow the universal law S <~ | log 2 £, 
where c is the central charge characterizing a correspond- 
ing universality class in conformal field theory (CFT) in 
the continuum limit. First obtained for continuous fields 
within the CFT framework Q , the result has by now been 
widely verified in discrete systems such as spin chains Q , 
harmonic chains [B[ , and fermions d, Q . 

That many-body systems should also manifest prop- 
erties of genuine multipartite entanglement (MPE) has 
been evidenced, for example, in spin chains [1(1 EH EH 
HH , and as shown in the simulation of many-body ground 
states within the matrix-product state framework 14) , 
MPE at criticality may be highly non-trivial. Connec- 
tions between critical behavior and MPE in ID spin mod 
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els have been particularly elaborated by 
which find non-analytic behavior of MPE at criticality 
and other signs of universality. 

In this Letter we further explore the MPE of critical ID 
systems by addressing a question in the spirit of renor- 
malization that naturally arises in this context: Given 
the scale-free nature of the continuum limit of the crit- 
ical system, how does the MPE between regions of the 
same size scale under coarse- or fine-graining, that is, as 
we vary the size of the regions (Fig. 1)? We investigate 
this problem in both harmonic and XY spin chains, using 
as a measure of MPE the geometric measure of Wei and 
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FIG. 1: Coarse graining of a system to N regions of size I, 
with Nl fixed. 



Goldbart Our main finding is that with respect 

to a partition of the system into equal regions of size £, 
the geometric MPE per region, £ , shows the logarithmic 
scaling behavior 
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at criticality. Here, c* is a constant that most probably 
is the central charge of the relevant universality class, as 
evidenced from our numerical results and the correspon- 
dence with the case bipartite case, for which a connection 
can be established independently using previous results 
in the literature. 

For an entangled state of N parties, the geometric en- 
tanglement measure [19| is defined as E(ip) — — log |A| 2 , 
where |A| 2 is the maximal overlap \(ip\4>i(f>2 ■ ■ ■ 4>n)\ 2 over 
all possible product states \^>i(f>2 ■ ■ ■ ^jv)j the optimal 
states and A are in turn solutions to the nonlinear eigen- 
value equation 
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(2) 



where <j>i stands for the exclusion of the state \4>i)- Here, 
we will be interested in the ground state \ip) of large 
circular chains of spins and oscillators, with the N par- 
ties representing regions of size £ each, so that the total 



2 



number of spins or oscillators is N£. Considerable simpli- 
fications of the nonlinear eigenvalue problem under these 
conditions allow for a feasible numerical implementation 
of the solution: First, due to translational invariance, 
the optimal state takes the form \<f>)® , and the prob- 
lem reduces to finding a single state \cf>) for one region. 
Furthermore, ground states of oscillator and spin chains 
in the XY model arc both dcscribable in the language of 
bosonic or fermionic gaussian states, where the connec- 
tion in the XY case is through a Jordan- Wigner trans- 
formation to effective fermion variables. Since partial 
tracing is a gaussian operation, a solution to Eq. ((2|) 
can be found within the class of gaussian states. As is 
well known, gaussian states can be characterized by a 
single covariance matrix[2(| that takes the block form 
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boson 



G K 
K T H 



M 



fermion — 



A C 
C T B 



(3) 

where Mboson is symmetric and Mf crm i on is antisymmet- 
ric. Further simplification follows from the fact that 
in the cases of interest, the conditions K = and 
A = B = are satisfied. For pure states, these con- 
ditions imply that in the bosonic case GH = 1, while in 
the fermionic case CC T = 1; thus, all the information of 
the gaussian state can be encoded in a single block (e.g., 
G or C) . We shall therefore refer to a single N£ x N£ 
matrix fi for the full state \tp) and an I x I matrix u; pt 
for the optimal state \(f>), where both i~2 and uj are either 
symmetric (bosonic) or orthogonal (fermionic). Finally, 
due to translational symmetry, has a Toeplitz form at 
the level of £ x I blocks, and may thus be brought through 
a matrix version of Bloch's theorem to the block-diagonal 
form uj = Q) r] u> v where the set {e" 1 } are the iVth roots 
of unity, and the uj v are hermitian (bosonic) or unitary 
(fermionic) £ x £ matrices, given by 



1 N ~ 1 

m=0 



(4) 



The optimal gaussian solution to the eigenvalue problem 
of Eq. [5] can then be recast as the equation 

1 o 

(5) 



°P N ^ oj op + u> n ' 

With this form, a numerical computation of u> pt can be 
obtained iteratively starting with a trial a; op t in the right 
hand side and iterating until a fixed point is reached. 
Our experience shows that in most cases twenty or so 
iterations suffice to reach an acceptable fixed point. 

Once u; opt is determined, the geometric entanglement 
can be obtained from the sum E(i(j) = E v , where 



£L(V0 = ±log 2 



det \ [u: op + uj n ] 



det oj op dct uj v 



(6) 



are the partial contributions from each sector and where 
the minus sign applies to the fermionic case. The av- 
erage of the partial entanglements E^tji) gives the en- 
tanglement per region or entanglement density E (tp) = 
E(ip)/N, which is our quantity of interest. 

We first consider the circular harmonic chain described 
by a Hamiltonian with a single parameter < a < 1, 

2 2 

H a = % + "9 — cSiSi+i) with i periodic in N£. Here, 
has entries given by = 2g(\i — j\) where g(l) is 

the momentum correlation function {popi) and is diago- 
nalized by discrete circular wave normal mode functions 
indexed by an angle 9 k — j^k, with eigenvalues given by 
the dispersion relation uj(6k) = y/l — a cos(9k). For fi- 
nite a, the correlation function decays exponentially with 
I with correlation length £ ~ f / y/2{l — a). Criticality 
corresponds to the limit a — > 1, in which the system be- 
comes gapless, (oj(9) oc \6\ for small theta), the correla- 
tion length diverges, and the correlation function exhibits 
power-law behavior g(l) ~ l/l 2 . 

It will be instructive to briefly develop a picture of the 
optimal solution based on the Bloch decomposition of O. 
In the harmonic chain, the oj^ are given in terms of the 
dispersion relation to (9) and circular plane waves with 
shifted node numbers 



k=0 



(7) 



where v„ — (j) mod 2tt)/2it is the winding number of rj. 
We can interpret oj^ as (twice) the momentum correlation 
matrix (7Ti7rj)^ for the vacuum state of a complex scalar 
field on a linear lattice of £ points, with Hamiltonian 
H n = ^tt' • 7r + ^(fy~V TI (j) and potential matrix = u>f , 



where 



i-j|>l 
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(this result is easily derived by applying the Bloch de- 
composition to the potential matrix for the whole chain). 
The case rj = corresponds to a translationally invari- 
ant closed chain of £ oscillators, while for 7ymod27r ^ 0, 
translational invariance is lost by the appearance of a 
"twisted" coupling between the first and last oscilla- 
tors, which in the continuum correspond to the twisted 
boundary conditions 0(0) = e tTI <j>(£) for a complex Klein- 
Gordon field. Now, due to periodicity in n, the uj v de- 
scribe a closed loop in the space of Hermitian £ x £ ma- 
trices, and uj op corresponds to a generalized "center of 
mass" for this loop with respect to the distance mea- 
sure (j6|). As £ becomes large, the perturbation to the 
free modes on the circular chain becomes noticeable only 
at the longest wavelengths, for which SX/X = vjk is ap- 
preciable, and gives rise to small ~ l/l corrections to the 
spectrum. Thus we expect that the distances E n should 
converge to a common value for large £. However, it is 
important to note that as a — * 1, the partial contribu- 
tion E v= q diverges independently of u; opt , this is due to 
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FIG. 2: Block size dependence of the geometric entanglement 
per block £ = E/N, as measured from its value at £ — 4 
(S£(£) = £(£) — £(4)) vs. log 2 I, for three values of the cor- 
relation length £ in the harmonic chain, for N = 10 blocks. 
Inset: Convergence of the slopes of partial contributions E v 
as function of log 2 £, for N = 20, at criticality. Curve or- 
der (starting from bottom curve) proportional to 1 77 — -zr | with 
T) € [0, 27r); dashed line corresponds to c/12 ~ 0.0833. 



the vanishing determinant of u) v =o and in turn to the 
dispersion relation u>{6) = for the 9 = mode at crit- 
icality. Thus, for 77 = 0, the relevant distance as a — > 1 
should in fact be taken to be E^Zq = E n= o — E^iv , where 
i?div = — ^logw(O) ~ ^og 2 £, is the divergent contribu- 
tion of the zero mode. The renormalized distances are 
well-defined at criticality and are indeed seen to converge 
to common values as £ — > 00. 

Turning now to our numerical results, we first note 
that the use of an entanglement density is indeed appro- 
priate, as the (renormalized) densities rapidly converge 
to iV-independent values at small values of N. This can 
be attributed to the fact that increasing the density of 
points in the previously mentioned closed loop of her- 
mitian matrices, has only minor effects on the location 
of the "center of mass". We also comment on the op- 
timum matrix u) op : the corresponding potential matrix 
V op = ijJ 2 ov was found to have, up to small fluctuations, 
the same structure of that of a chain of length £ with the 
same value of a but with open boundary conditions; thus, 
in a first approximation, the optimum potential matrix 
appears to be the average of the matrices as one may 
expect from the center of mass picture. Turning then to 
the entanglement density as a function of £, in Fig. 2 we 
present our results for different correlation lengths £ . As 
can be seen, the increase in entanglement for non-critical 
values saturates when £ ~ £. This limiting behavior can 
be attributed to the entanglement of oscillators within 
a distance ~ £ from the interfaces between regions, in 
which case the total entanglement is determined only by 
the number of partitions. On the other hand, at critical- 
ity £ = 00 1 the- renormalized entanglement per region is 
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FIG. 3: Block size dependence of the geometric entanglement 
per block £ = E/N, as measured from its value at £ = 4 
(S£(£) ee £ (£)-£ (4)) vs. log 2 £, for three values of (7, A) on the 
phase plane in the XY model, for N — 20 blocks. From top 
to bottom: vicinity of the bosonic critical line (XX model), 
fermionic critical line (Ising) , and a non critical value showing 
saturation. Inset: Convergence of the partial contributions 
E v as function of log 2 £, for N = 10, at the Ising model critical 
point. Curve heights follow same order as in Fig; dashed line 
corresponds to c/12 ~ 0.0417 . 

found to scale with £ as 

£<™»)(V0 = K*log a * (9) 

where k* is a coefficient that one may expect to depend 
on the CFT central charge c = 1, as does the bipar- 
tite entanglement entropy. For the finite values consid- 
ered in our computations, the coefficient k* in fact varies 
slowly as £ is increased and it therefore becomes diffi- 
cult to infer the limiting value from the total entangle- 
ment density. However, an examination of the instant 
slopes of the partial contributions E n reveals (as shown 
in Fig. ) that they all converge to a single limiting value, 
from both directions. We have estimated this limiting 
value by fitting the instantaneous slopes to the function 
K>u{£) = k* + A(y)l~ x , where A(y) is quadratic in v and 
symmetric about v = 1/2 (from periodicity), for vari- 
ous values of partitions up to N = 100, and obtain that 
k* = 0.0837 with a 5% error, independently of the num- 
ber N of partitions. Within the margin of error, this 
value is consistent with k* = c/12 where c is the central 
charge c = 1 of the bosonic CFT universality class. 

To test the universality of our results, we consider 
the XY spin chain model, which includes two criti- 
cal regions associated with different CFT universality 
classes. The XY Hamiltonian takes the form H\ n = 
- £ XZ, + ^X t X l+1 + ^YiY i+1 , where Xi,Y it z\ are 
Pauli matrices at the site i, and the anisotropy 7 and 
the field strength A parametrize the ground state phase 
diagram. The critical lines are 7 = 0, associated with 
the bosonic (c = 1) universality class, and A = 1, as- 
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sociated with the fermionic (c = 1/2) class. By means 
of a Jordan- Wigner transformation, H\ n can be trans- 
formed to a free fermions system, the ground state of 
which is gaussian fermionic and can be characterized by 
an orthogonal matrix fi with eigenvalues 



u(9) 



(cos(0) - A) -i7sin(0) 
v/(cos(#) - A) 2 + 7 2 sin(0) 



(10) 



and circular plane waves with momentum label 8 as 
eigenstates. The corresponding uj v matrices can be ob- 
tained by using this expression for uj(9) in Equation [3 

As depicted in Figure 3, the XY model also shows sat- 
uration for non- critical regions and logarithmic scaling 
behavior for spin chain models in at criticality. For the 
critical Ising point, 7 = A = 1, very good agreement is ob- 
tained with k* ~ 0.043, consistent with c/12 for c = 1/2. 
In the critical XX case, A = 0, < 7 < 1, our itera- 
tion scheme converges slower, and large £ values become 
harder to obtain numerically due to precision losses at 
every iteration. Still, logarithmic scaling was evidenced 
with k* consistent with c/12 for c — 1 within 15%. 

The slope value c/12 is consistent with the results for 
N = 2, which are analytically tractable, owing to the 
fact that in the bipartite case the geometric entangle- 
ment is the logarithm of the largest Schmidt coefficient 



of the state. As discussed in (2ll . [22J, with respect to 
any bi-partite split the ground state of a quadratic bo- 
son or fermion hamiltonian can always be expressed as 
the product of two-mode entangled states of the form 
(1 ± e~ l3 ) ±1 l 2 ^2 n e~ (3n \n,n) where the /3's are related 
to the symplectic eigenvalues of the reduced covariance 
matrix and n goes from zero to 00 for bosons, or 1 for 
fermions. The geometric entanglement is therefore the 
"free energy" E = ±^olog(l ± e _/3 ) coming from the 
n = coefficients. In the continuum limit this quantity 
is equal to half the von Neumann entropy, if the den- 
sity of modes is constant for the range of contributing 
modes /3 < 1, as^curs at criticality (a rigorous deriva- 
tion is found in 
Cardy 



23])- Using the result of Calabrese and 



24j . where for a symmetric split of a chain of 



size 2£, the entanglement entropy is found to scale like 
S ~ I log £, the resulting geometric entanglement density 
is then found to be £n=2 = hE = jS, and thus to scale 



as 



12 



logi 



Elsewhere, Bravyi 25] has proposed a generalization of 
the entanglement entropy for multipartite states, corre- 
sponding to the minimum attainable Shannon entropy of 
the joint measurement outcome distribution when con- 
sidering all possible multilocal measurement bases. Since 
the overlap between an entangled state and any prod- 
uct state cannot exceed in magnitude the maximal value 
|A| from eq. ([2]), the Shannon entropy of the squares 
of the coefficients in any product basis expansion of the 
state cannot exceed the entropy — log 2 | A| 2 of an equally- 
weighted superposition of 1/|A| 2 product-state terms. 
The geometric entropy is therefore a lower bound on the 



multipartite entanglement entropy. Roughly speaking, 
such an entropy may be considered a measure of the ef- 
fective Schmidt number [3], or number of terms in a 
multilocal orthogonal decomposition of the state. Our 
result therefore suggests that at criticality this effective 
number must scale with the region size I no slower than 

Nc . . . 

~ 1 12 . An interesting open question is the extent to 
which the exponent ^ is a good characterization of the 
actual effective Schmidt number of ID ground states for 
critical systems. 

To conclude, in the present work we found convincing 
evidence that multi-partite entanglement in ID systems, 
manifests at criticality a logarithmic scaling behavior as 
well as universality; properties which have so far been es- 
tablished only for bi-partite entanglement at criticality. 
This finding, has been here demonstrated for free Gaus- 
sian harmonic and spin chain models, and will hopefully 
be further tested in other discrete solvable models, as 
well as in the framework of conformal field theory. It 
would be also very interesting to formulate the problem 
at hand in conformal field theory, since the connection 
of the present result with some properties of CFT is at 
present an intriguing open question. Finally we hope 
that, the emerging understanding of the properties and 
structure of entanglement in many-body systems, will 
also turn helpful in further developing tools for study- 
ing many-body systems. 
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